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ABSTRACT 

We present a simple method, based on the deformation of spherically sym¬ 
metric potentials, to construct explicit axisymmetric and triaxial MOND density- 
potential pairs. General guidelines to the choice of suitable deformations, so that 
the resulting density distribution is nowhere negative, are presented. This flexible 
method offers for the hrst time the possibility to study the MOND gravitational 
held for sufficiently general and realistic density distributions without resorting 
to sophisticated numerical codes. The technique is illustrated by constructing the 
MOND density-potential pair for a triaxial galaxy model that, in the absence of 
deformation, reduces to the Hernquist sphere. Such analytical solutions are also 
relevant to test and validate numerical codes. Here we present a new numerical 
potential solver designed to solve the MOND held equation for arbitrary den¬ 
sity distributions: the code is tested with excellent results against the analytic 
MOND triaxial Hernquist model and the MOND razor-thin Kuzmin disk, and a 
simple application is hnally presented. 

Subject headings: gravitation — stellar dynamics — galaxies: structure — meth¬ 
ods: analytical — methods: numerical 
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1. Introduction 


Milgrom (1983) proposed that the failure of galactic rotation curves to decline in Kep- 
lerian fashion outside the galaxies’ luminous body arises not because galaxies are embedded 
in massive dark halos, but because Newton’s law of gravity has to be modihed for helds 
that generate accelerations smaller than some characteristic value Oq. Subsequently, in order 
to solve basic problems presented by this phenomenological formulation of the theory (now 
known as Modihed Newtonian Dynamics or MOND), such as conservation of linear momen¬ 
tum (Felten 1984), Bekenstein & Milgrom (1984) substituted the heuristic 1983 model with 
the MOND non-relativistic held equation 


V- 


niv^ii\ 


dv 


Uo / 


V0 


AttGp, 


( 1 ) 


where ||...|| is the standard Euclidean norm, 0 is the gravitational potential produced by the 
density distribution p, and V0 —0 for ||x|| —cx). As stressed by Bekenstein & Milgrom 
(1984), the equation above is obtained from a variational principle applied to a Lagrangian 
with all the required symmetries, so the standard conservation laws are obeyed. Thus, 
equation (1) plays in MOND the same role as the Poisson equation 


V^0N = 4:7rGp 


( 2 ) 


in Newtonian gravity, and the MOND gravitational held g experienced by a test particle is 


g = -V0. 


( 3 ) 


As well known, in the regime of intermediate accelerations the function p is not fully 
constrained by theory or observations, while in the asymptotic regimes 


pit) ~ 


t for t 1, 
1 for t 3> 1. 


( 4 ) 


Throughout the paper we conform to the standard assumption 


pit) 


t 

TtW 


(see however Famaey & Binney 2005). 


( 5 ) 


From equation (4) it follows that equation (1) reduces to the Poisson equation when 
||V0|| 3> Oo, while the limit equation 


V-(||V0||V0)=47rGaop, 


( 6 ) 
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obtained by assuming yu(t) = t in equation (1), describes systems for which (or regions of 
space where) || V0|| Oq, i.e. systems for which the MOND predictions differ most from the 
Newtonian ones. As a consequence equation (6), characterizing the so-called ‘deep MOND 
regime’ (hereafter dMOND), is of particular relevance in MOND investigations. 

Nowadays a considerable body of observational data seems to support MOND well 
beyond its originally intended held of application (see, e.g., Milgrom 2002; Sanders & 
McGaugh 2002), making this theory an interesting alternative to the Cold Dark Matter 
paradigm. It is thus natural to study in detail MOND predictions, in particular focusing on 
dMOND systems, i.e. systems that should be dark matter dominated if Newtonian gravity 
holds. Potential problems of the theory have already been pointed out by various authors 
(see, e.g. The & White 1988; Buote et ah 2002; Sanders 2003; Ciotti & Binney 2004; Knebe 
& Gibson 2004; Zhao et ah 2005), but further analyses are needed to reach hrmer con¬ 
clusions. Unfortunately, MOND investigations have been considerably slowed down by the 
lack of aspherical density-potential pairs to test theory predictions in cases more realistic 
than those described by spherical symmetry: the search for a sufficiently general method to 
construct aspherical MOND solutions is the subject of this paper. 

In MOND, the main difficulty to obtain exact aspherical density-potential pairs (or to 
build robust numerical solvers) originates from the non-linear nature of the theory, which 
makes impossible a straightforward use of the analytical and numerical techniques available 
for the Poisson equation (such as integral transforms or expansion in orthogonal functions). 
In addition, a simple relation between the Newtonian and the MOND gravity helds in gen¬ 
eral does not exist. Although equation (2) can be used to lower the order of equation (1) by 
eliminating the source density, it follows that /i(|| V0||/ao)V0 and V0 n differ by some un¬ 
known solenoidal field^ S = curlh. Remarkably, Brada & Milgrom (1995; hereafter BM95) 
showed that when the modulus of the Newtonian gravitational held produced by p is a 
function of 0 n only, S vanishes, and so the MOND acceleration g is related to gN = —V0n 
by 

h( —)g = gN, (7) 

V«o/ 

where g = ||g||. Equation (7) coincides with the original MOND formulation of Milgrom 


^That the solenoidal field S in general cannot be arbitrarily set to zero is due to the fact that the MOND 
acceleration field must be derived from a potential, while equation (7), if correct in general, would imply 
that V||g|j A g = 0, an identity which is not necessarily true when g is derived from a potential. 



(1983) and can be solved algebraically: particularly simple cases described by this equa¬ 
tions are those in which the density distribution is spherically or cylindrically symmetric, 
or stratihed on homogeneous planes. In such cases the MOND potential is not more diffi¬ 
cult to construct than the corresponding Newtonian potential. This is the reason why all 
the (astrophysically relevant) analytical MOND density-potential pairs known are spheri¬ 
cally symmetric: in fact, the only exact, aspherical MOND density-potential pair available 
is the axisymmetric razor-thin Kuzmin disk (and the derived family; BM95), for which 
S'N = 5'n(0n)- In all the other cases one is forced to solve numerically equation (1). 

The importance of a general method to obtain the explicit MOND potential of den¬ 
sity distributions with prescribed shape and stratihcation is then obvious, because it would 
allow for orbit integration, without resorting to numerical integration of equation (1); in 
addition, analytical solutions could be used to test numerical MOND solvers in more real¬ 
istic cases than spherical symmetry. In this paper we show that such an approach can be 
devised. In particular, we show how a “seed” spherical density distribution can be deformed 
in an axisymmetric or triaxial density distribution with analytical potential satisfying the 
MOND equation, by means of a pair of very simple existence theorems and (for example) 
by using building blocks obtained from the Ciotti & Bertin (2005; hereafter CB05) method. 
In addition, as a complement to the analytical method, we also illustrate a new numerical 
MOND solver based on spectral methods, which adds to the short list of the others available 
(Milgrom 1986; BM95; see also Brada & Milgrom 1999). 

The paper is organized as follows. In Section 2 we present the general method (post¬ 
poning to the Appendix the proof of three technical results on which the method is based), 
and in Section 3 we construct families of triaxial prohles obtained by deformation of the 
Hernquist (1990) sphere. The original numerical code developed to compute the MOND 
potential of generic density distributions is then described in Section 4, where we also show 
how well it recovers the analytical triaxial potentials of Section 3 and the MOND Kuzmin 
density-potential pair. The code is hnally used to investigate the field S of highly flattened 
triaxial density distributions. The main results of the paper and possible future applications 
are hnally summarized in Section 5. 
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2. The general method 


Two different approaches are possible to construct explicit solutions of equation (1); in 
the most obvious one attempts to recover the potential of a given density p. Unfortunately, 
no explicit aspherical density-potential pairs (with the exception of the Kuzmin disk) have 
been obtained so far following this approach. In addition, we note that even “innocent” 
density distributions can produce puzzling behaviors of the MOND acceleration field at 
special places, as illustrated by the perturbative case analyzed by Ciotti & Binney (2004). 

In the alternative “0-to-p” approach - studied in this paper - one determines the den¬ 
sity held by application of the MOND operator to a prescribed potential. This approach 
is straightforward because only differentiation is involved, but negative densities can be 
produced if the potential is chosen without the necessary care, as well known also in the 
Newtonian case. Thus, a hrst problem posed by this approach is to hnd a criterion to choose 
0 such that p is positive. In addition to the positivity condition, one would like to be able 
to control the shape and the radial trend of the resulting density distribution: for example, 
in Newtonian gravity homeoidally stratihed potentials lead to curious toroidal density dis¬ 
tributions (e.g. Binney 1981, Evans 1994). As we will see, it is possible - to some degree - 
to answer positively to these issues. 


The idea behind our method is to compute the MOND potential of a spherical density 
distribution. The obtained potential is then mapped on a new spherical density by means 
of the Laplace operator. This new density is then suitably deformed and the associated 
(Newtonian) potential is determined by standard methods. Finally, the MOND differential 
operator is applied to the deformed potential, and the hnal density distribution is obtained. 
In practice, the method consists of the following steps. 

1) We start by selecting a spherical density distribution po{r) with the desired radial be¬ 
havior: for example, it could be a y-model (Dehnen 1993; Tremaine et ah 1994), a King (1972) 
density distribution, or a classical polytrope such as the Plummer (1911) sphere. We then 
calculate its dMOND potential (l>o{r) by using equation (7) in spherical symmetry: 

d(j)o A/GMo(rj^ 

-r = 9o = -, ( 8 ) 

dr r 

where Mo{r) = An po(r)r^dr. As is well known, </)o(r — > cx)) ~ V GM^aQ In r for any mass 


first attempt to construct the MOND potential of an oblate galaxy along this line was carried out by 


Hongsheng Zhao (private communication). 
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distribution of finite total mass Mq. We stress that we are calculating the dMOND potential 
of po irrespective of whether the density distribution is actually dMOND everywhere: for 
example, the Jaffe (1983) model near the center is not dMOND because its acceleration held 
diverges there, yet we can integrate equation (8) to obtain 0o- The reasons for doing so are 
that equation (8) is much easier to integrate than the corresponding MOND (spherically 
symmetric) equation, and that Theorem 2 assures that both the Laplacian and the MOND 
operator applied to 0o will produce positive densities. In particular, the Laplace operator 
applied to 0o leads to the density distribution 


Pno 



pQ{r)r ^ -jMoir) 
2^JMo{r) 47rr2 


(9) 


while the MOND operator in equation (1) produces 


Pmo 



Utt- I 

V«0/ 

r° + 47rG(l + PoV«g)J 


( 10 ) 


By construction 0o is the dMOND potential of po, the Newtonian potential of pno, and the 
MOND potential of pmo- Obviously, pmo coincides with po where the model is in dMOND 
regime, and with pNo where the gravity held go is strong enough. From a more quantitative 
point of view, Pno(^ oo) oc when the total mass of po is hnite, and Pno(^ —0) oc 
7 >-(i +“)/2 jf oc (a < 3) in the central regions. This is apparent from Fig. 1, where we 
plot Po (solid lines), pno (dashed lines) and pmo (empty symbols) for 7 -models with 7 = 0, 
1, and 2. 


2) In the second step of the procedure we look for an aspherical density pni to be added 
to Pno so that pno + ^Pni is positive for some A > 0, and (dehned by V^0i = dvrCpNi) 
is explicitly known. The associated total (Newtonian) potential is thus 0 = 00 + A0i. The 
main reason for working with Newtonian potential at this stage is that the properties of 
density-potential pairs obeying the Poisson equation are well known, and so we can work 
with conhdence and control the deformations of the total pno + ApNi. 

3) The hnal step is the explicit evaluation of the dMOND operator in equation (6) for 
the potential 0 = 0o + A0i, and this can be done with computer algebra packages such as 
Maple or Mathematica. For A = 0 one recovers 0o (which produces the positive po), and for 
continuity one could expect that a positive p (with a radial trend similar to that of po) is 
also obtained for sufficiently small A, though we will show that this is not always the case^. 


^See the discussion below equation (15). 
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log r/r^ log r/r^ log r/r^ 


Fig. 1.— The Newtonian (dashed lines) density prohle pno [equation (9)] and the full MOND 
(empty symbols) density prohle pmo [equation (10)] derived from the dMOND potential of 
the spherical y-model Po = (3 — 7 )Morc/[ 47 rr'''(r + rcY~^] for 7 = 0,1, 2 (solid lines). In all 
cases we adopted GMq = lOOaor^. 
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In general, for given p and pNi a critical value of A exists such that only smaller values A 
correspond to physically acceptable (i.e. nowhere negative) p. If the resulting density is 
positive, then Theorem 1 assures us that also the full MOND operator applied to 0 will 
produce a positive density distribution (such that its limit for A —0 is Pmo), provided that 
Pno +ApNi > 0. Note that the last requirement is just a sufficient condition for the positivity 
of p. 

A hrst hint about the positivity and the shape of the resulting density can be obtained 
by studying the linearization of the dMOND operator. In fact, it is easy to prove that the 
A-linearized density associated with 0 = 0o + A0i is 

4^ 90 + (4xG9n„ + + 0(A=), (11) 

where the spherical symmetry of 0 o has been exploited by using the standard spherical 
coordinates (r, d, <p), and it is intended that = (f)i{rffi,ip) and go ^ 0 . 


From the discussion above it should be obvious that the delicate step in the procedure 
is the choice of pni such that pno + ApNi is nowhere negative, 01 is analytical, and the 
dMOND operator produces the sought deformation on p. Simple choices satisfying the hrst 
two requests could be the addition of Miyamoto & Nagai (1975) or Satho (1980) disks: 
however, an even more general approach can be devised. In fact, families of explicit and 
easy-to-calculate Newtonian (p^ii^i) pahs can be obtained from the homeoidal expansion 
technique (e.g., see CB05). In practice, for the present problem one can adopt as pNi the 
linear term in the homeoidal expansion of a seed density g, i.e. 

Pm = {ey^+ Vz^)^-^, (12) 

r 


where the dimensionless parameters e and p are the battenings of the expanded homeoidal 
density distribution. The potential 0i generated by the distribution pNi is written in terms 
of one-dimensional radial integrals that usually can be evaluated explicitly (see equations [5] 
and [ 6 ] in CB05): 

-p- = (e + p) 0 i(r) + (ep 2 + pz2)02(r), ( 13 ) 


where 


01(r) = 


Q{m)w? 


^ nr\ ^ 2 

^-3;5r’"+3 


g{m)mdm, 


02(r) = 


g{m)rffidm. 
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We stress that in the present application e and rj are just linear parameters, so without 
loss of generality the multiplicative coefficient A can be considered to be contained in them, 
and their value is not restricted by the limitations described in appendix A of CB05. Note 
also that pNi in equation (12) is negative when it is derived from an oblate axisymmetric (the 
e = 0 case) or triaxial g, and so particular care is needed in the choice of g, in order to have 
Pno + Pni > 0. However, for a g that can be approximated by a power-law in its external 
regions, pni(^" —^ cxd) ~ —g, and so any hnite-mass g produces through equation (12) a 
positive total density at large radii, where pno ~ Finally we remark that, for pNi as in 
equation (12), the quantities d(pi/dr and d'^^’i/dr'^ in equation (11) are particularly simple 
because, according to equation (13), the angular part of is just a multiplicative factor of 
the radial function r^'^ 2 - 

We note that the seemingly obvious choice of taking pNi to be the homeoidal expansion 
term derived from does not work. This is because pno ~ at large radii, so pni oc 
— (esin^ dsin^ (p + t] cos^ d)/r‘^ and 01 oc —2(e -|- p) lnr/3 -|- (e sin^ d sin^ <P + P cos^ 'd)/3 for 
r —oo: equation (11) then reveals that the density is negative for sufficiently large r in the 
conical region cos P < 1/v^. As we will see, the regions near the 2 : axis are in fact the most 
delicate when applying the analytical technique presented in this paper: unphysical models 
usually develop negative densities near the 2 :-axis. 

We conclude this Section by noting that, when truncated at the first order in A, 0o + A0i 
and the A-linearized p in equation (11) allow for the explicit evaluation of the A-linearized 
solenoidal held S = ||V0||V0/ao — V0 n = ASi -|- O(A^) (the zeroth-order term in A is 
obviously missing because it refers to the spherical density component). If 0i is given by 
equation (13) then p in equation (11) belongs to the family considered in CB05 and its 
Newtonian potential can be derived in closed form (while in general this is not possible 
when considering the density obtained by the application of the dMOND operator without 
A-linearization). 


3. A simple application: triaxial Hernquist galcixy 

In order to illustrate the method described in Section 2 we now construct a triaxial 
galaxy model with a realistic density distribution and analytical MOND potential. 

Following Step 1, we start with the spherically symmetric Hernquist (1990) density 



Fig. 2.— Isodensity contours (left panels) and density profiles (right panels) for two ana¬ 
lytical dMOND axisymmetric (e = 0) Hernquist models with 01 as in equation (21). The 
density prohles are taken along a radius in the equatorial plane (dashed lines) and along the 
symmetry axis (solid lines). The model in the top panels has 0 = 5 and fj = 0.01, while 
the model in the bottom panels has j3 = 5 and fj = 0.02. 
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Fig. 4.— The dMOND solenoidal field S = gg/ao — gN (normalized to in the meridional 
plane for three different axisymmetric Hernquist models (left panels). In the right panels 
lines of constant HSH/^fN are shown. From top to bottom; the model of Fig. 2 with fj = 0.02, 
the model of Fig. 3 with g = 0.4, and the homeoidally stratified axisymmetric Hernquist 
model [equation (37)] with e = 0 and g = 0.5. 
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distribution 


Po{r) = 


Ma 


(16) 


2nr^ s(l + s)^ ’ 

where s = r/Vc, we recall that this density distribution, when projected, is a good approxi¬ 
mation of the law (de Vaucouleurs 1948) over a large radial interval, thus providing a 
realistic description of the luminous density distribution of elliptical galaxies. From equa¬ 
tion (8) go = y/GMotto/rc{l + s), and so the dMOND potential is 

00 = GMotto ln(l -|- s), (17) 

while from equation (9) 


(18) 

r~^ in the central regions and 


y/GMotto s -\- 2 
AnGrl s(s-|-l)2' 

In agreement with the discussion in Section 2, pNo ~ Po 
Pno oc in the outer parts of the dMOND Hernquist model . 

According to Step 2, we look for a density distribution p^i to be added to pno such that 
the total density is positive and the potential 0i is explicitly known. A family of densities 
that can be easily expanded with the CB05 technique and admit a simple analytical potential 
01 is 


p(r) = 




[s + PY 


(a > 3), 


(19) 


where go is a density scale, and 0 is the scale-length in units of Vc- There is nothing special 
about the distribution (19), which has been chosen only by simplicity arguments: note 
however that at both small and large radii the total density is dominated by pNo, so its 
positivity has to be checked only for the intermediate regions. Here we restrict to the case 
a = 6 (see also Muccione & Ciotti 2004), and from equation (12) 

g (e sin^ d sin^ + V cos^ P)s 


Pm = -6po0 


{s + PY 


so that 


01 = AirGgoP^r, 


5^2 

C 


(e -I- 7]) (s^ -I- 30s -f 0^) ^ (e sin^ d sin^ <p + p cos^ d)s^ 


3O02(s + 0)3 


5(s -F 0)^ 


( 20 ) 


( 21 ) 


Before computing the dMOND operator for 0o + 0i it is worth studying its linearization 
for e —0 and p —0, and from equation (11) we obtain 


+ P%e + p) 


prl 1 

Mo 27rs{s + 1)^ 

60® (e sin^ d sin^ <p + p cos^ d) 


s® + 7ps^ + 902s® + 02(0 - 6 )s 2 - 20®(0 + 6)s - 60^ 
15 (s + 1)2(s + 0)7 

6s® - (230 - 3 )s2 + 0(0 - 24)s + 302 
15(s+l)2(s + 0)7 ’ 


+ 


( 22 ) 
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where the new dimensionless parameters e and fj are dehned as e/e = fj/rj = gor/^/G/M^a/. 
Thus, p ~ po ~ for r —0, while for r ^ oo 


prl 

Mo 



e + T] 
15 


1 , „ ^5 e sin^ d sin^ <P + P cos^ d 


^ + 12 / 9 ^ 


5s6 


(23) 


i.e., the linearized density retains the spherical symmetry and radial trend of po for r —0 
and r —cx), so negative values of p can be present at intermediate distances only from the 
center. For instance, restricting to the axisymmetric case (e = 0) with j3 = 5, according to 
equation (22)negative values of p first appear for fj > 0.018 along the 2 :-axis near 2 : ~ 4rc, 
and smaller values of fj correspond to an everywhere positive p. 


The formula resulting from the evaluation of the dMOND operator is not reported 
here, and we limit ourselves to show the isodensity contours in a pair of representative 
axisymmetric (e = 0) cases, hxing j3 = 5. In Fig. 2 (top left panel) we show the isodensity 
contours in the meridional plane for the model with fj = 0.01. Note how the density is 
spherically symmetric both in the central regions and at large distances, in accordance with 
equation (22). It is particularly important to stress that the maximum flattening of the 
density corresponds to an axis ratio ~ 0.8, even though fj is much smaller than 0.2. This 
at variance with the Newtonian case (in which the flattening of the expanded density is the 
same as that of the seed density for small deformations), and this is a result of the non-linear 
nature of MOND. In the top right panel of Fig. 2 we show the density prohles along the 
symmetry axis (solid line) and along a radius in the equatorial plane (dashed line) for the 
same model. It is apparent how the radial trend of the spherical Hernquist model has been 
nicely preserved by the imposed deformation, again in accordance with equation (22). Note, 
however, how the density prohle along the axis shows a small depression, which becomes 
more apparent in the model fj = 0.02 (Fig. 2, lower panels). This last model is near the 
consistency limit in fj, and the axis ratio of the flattest isodensity surface is ~ 0.6. Larger 
values of fj would produce a negative density along the 2 : axis at 2 : ~ dr^, in agreement with 
the results of the linear analysis. 


It should be clear that the presented model is just one of the many analytical triaxial 
Hernquist models one can construct. For example, one could adopt the 0i derived from 
the homeoidal expansion of a power-law g (CB05), or just add a quadrupole potential [i.e., 
assume = 0 in equation (13)]. For instance, in Fig. 3 we show the density obtained by eval¬ 
uating the dMOND density of the potential 00 + 0i, where 0o is still given by equation (17), 
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1 77777—(esin^'(9sin^(^ + ?7 cos^-(9)s 

01 = ^/GM^ao- - . I ' -^ (24) 

[s + (iy 

The two models shown are characterized by e = 0, 0 = 1, ?7 = 0.2 (top) and r] = 0.4 
(bottom). Overall, we can note that for these models the isodensity surfaces are flat also 
in the central regions and in the models outskirts, and that the symmetry axis is again the 
region where the model becomes unphysical for too large rj. 

We conclude this Section discussing briefly the global properties of the dMOND held 
S = gg/do — gN for the presented models. In Fig. 4 (left) we show the held S/^^n in the 
meridional plane, while in the right panels we plot lines of constant S/g^, where S = ||S|| 
(The Newtonian held gN has been evaluated numerically with a spectral Poisson solver; see 
Section 4). The quantity S/g-^ gives direct information about the relative importance of the 
solenoidal held with respect to the Newtonian gravitational held. Remarkably, as already 
pointed out by BM95, also in our case the solenoidal held is typically almost everywhere much 
smaller than g-^: we hnd S'/^'n ~ 0.05 for the family of axisymmetric and triaxial models 
obtained from equation (21). However, we also hnd that this is not necessarily the rule: for 
example, S'/^'n ~ 0.25 in the central regions of the models obtained from equation (24), and 
other cases of non-negligible solenoidal held will be discussed in the next Section. 


4. The numerical MOND potential solver 


The results presented in the previous Sections, albeit encouraging, suggest that in order 
to have a full control on the density held producing the MOND potential the best tool is 
still the numerical solution of equation (1). To our knowledge the only MOND numerical 
solvers presently available are the two-dimensional code developed by Milgrom (1986), and 
the BM95 three-dimensional, multi-grid solver [used by Brada & Milgrom (1999) in their 
N-body code]. Here, as a complement to the analytical method, we present a new numerical 
three-dimensional MOND potential solver based on a spherical coordinates grid, which can 
be easily implemented in a standard particle-mesh N-body code. 


The goal is to solve numerically equation (1) for 0(x), i.e. the non-linear elliptic bound¬ 
ary value problem dehned by 


M[0(x)] = V 


/i 


V0(. 

do J 


47rGp(x) = 0, g = 0{r ) for r — oo, (25) 


where x = {x,y,z), r = ||x||, p(x) is assigned, and g = ||g|| with g = —V0(x). 
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The numerical scheme is based on the iterative Newton method, which is robust and 
relatively simple to implement. At the (n + 1) — th iteration we indicate with the 

increment function to be determined so that +50*^"^ is a better approximation of 

the solution, while and are quantities known from the previous iteration. 

In our case, is the solution of the linear equation 

= -M , (26) 

where the linear operator 

SMin) ^ y . ^ ^'(AgW ^g(n) . y^ j ^27) 


satisfies the identity 


M - M + O , (28) 


where /ao) and = p' /9^"'^a-o are known. Equation (26) is solved 

by inversion of and the procedure is repeated until numerical convergence to the 

solution 0num (dehued by max |M[0num]| < where e is a prescribed tolerance) is attained. 

The solution of equation (26) requires in general a hnite difference approximation of the 
spatial derivatives and the inversion of a three-dimensional matrix. Note that boundedness 
of the inverse of the operator assures quadratic convergence of the scheme for 0^°^ 

sufficiently close to the sought solution (e.g. Stoer & Bulirsch 1980). In our implementa¬ 
tion, designed for hnite mass distributions, we represent equation (25) and (26) in spherical 
coordinates, which allow for a simple assignment of the boundary conditions. Moreover, we 
approximate the linear operator (27) with 


= 


A 

dr 


{Ls + L, 


(29) 


where 

f =_LAT dA f =J_21 

sin -d (9-d \ dd) ' sin d dip'^ 


(30) 


are the angular components of the Laplace operator, = (l/47r) J sinddddip, 

and the variation o(; jg neglected. It is easy to show that with this new 

operator quadratic convergence in equation (26) is replaced by linear convergence since 


M - M + O [5(1)^^^' . 


(31) 
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However, the inversion of is much simpler because we are now in the position to use 

the full power of spherical harmonics. In fact, after expanding the source term M and 

the unknown function <50^”^ in spherical harmonics 


l,m 


(e.g. Binney & Tremaine 1987), the modihed equation (26) reduces to 


•9 f 2^(n) <9 

or \ or 


fl^^\r)l{l + 1) 


6(j)' 


(n), 

L,m 


r = 


-M 


l,m 


(33) 


involving derivatives only in the radial coordinates. The solutions S(j)^P^{r) are back trans¬ 
formed into 5(l)^'^\r,^,ip) according to equation (32), and the new components of the accel¬ 
eration held = g("-) _ are evaluated. 


In order to solve equation (33) we introduce the invertible mapping for the radial coor¬ 
dinate 

r(0 = tan"C, (0 < ^ < n/2), (34) 


so that the inhnite radial range is mapped onto the hnite interval [0,7r/2), 

d cos^ ^ d 

dr atan"”^ (9.^’ 


(35) 


and the boundary conditions for the radial derivative of the potential at .^ = 7 r /2 are au¬ 
tomatically satished (Londrillo & Messina 1990); in our applications we adopt a = 1 or 
CK = 2. In practice, the coordinates are discretized on a uniform grid 

of X X N^p points (7V^ = 64, = 32, and = 64, in typical applications), the d 

and derivatives are evaluated using the spectral representation of hnite order Legendre- 
Fourier polynomials, and the radial derivatives are approximated by second order centered 
diherences in the coordinate. Thus, for each (/,m) the discretized operator to be 

inverted is represented by a pentadiagonal matrix over We note that for / > 0 we solve 
equation (33) for with the boundary condition = 'k/2) = 0, while for / = 0 

equation (33) can be solved directly for the radial component of the acceleration increment. 


4.1. The tests 

We tested the numerical code using the analytical axisymmetric and triaxial MOND 
models described in Section 3. We assigned the analytical density distribution obtained for 
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Fig. 5.— R.m.s. (solid line) and maximum (dotted line) relative error on the acceler¬ 
ation (calculated over all the space) as a function of the number of radial grid points 
= 2N^ for two triaxial Hernquist models. Empty circles refer to the model with (pi 
in equation (21) with P = 5, fj = 0.02 and e = 0.01, while solid squares to the model with 
01 in equation (24) with P = 1, r] = 0.4 and e = 0.2. 

different values of aor"^/{GM q), e, t], and 0, and we recovered the corresponding numerical 
MOND potential, comparing it with its analytical expression. The accuracy of the numeri¬ 
cally estimated field gnum = — V0num is quantihed by the relative error ||g — gnum||/5', and 
in all the studied cases we found very good agreement between the analytical and numerical 
acceleration held, even for models near the consistency limit such as those in the bottom 
panels of Fig. 2 and 3. For example, in Fig. 5 we plot the r.m.s. (solid line) and maximum 
(dotted line) relative error calculated over all grid points, as a function of = 2N^ for 

two triaxial Hernquist models: one with (pi as in equation (21) (empty symbols), the other 
with (pi as in equation (24) (solid symbols). It is apparent how the typical errors decrease 
down to values ~ 10“^; such values are reached in 15 — 20 iterations, taking few seconds on 
a common PC. 

As a second set of tests, we studied the razor-thin Kuzmin disk (BM95). This is a quite 
severe test, because of the singular nature of the density distribution along the 2 : direction. 
In this case, we found numerical convergence to the solution, with relative r.m.s. error 
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5.6 X 10-3 for N^ = N^ = 128, and = 64. 

We note that onr code converges also when dealing with mnlti-centered density distri- 
bntions, for which the spherical grid is clearly not optimal. In snch cases the convergence is 
not as fast as in the case of density distribntions centered in the origin. For example, the 
MOND helds of gronps of 2, 3, and 4 Hernqnist spheres are compnted in 20 — 25 iterations, 
with tolerance e ~ lO”®, = 64, = 32, and N,p = 64. 


As an illustrative application of the code, we evaluated the dMOND held for the family 
of homeoidally stratihed triaxial Hernqnist models 


Mn 


P = 


2nr^{l — e)(l — ? 7 )m(l + m)^’ 


where 


9 tjU 

m = ^ + 


y 


+ 


r 2 (l — e)^ r^(l — 77 )^’ 


(36) 


(37) 


and we investigated their solenoidal held for diherent battenings. For example, the model 
in Fig. 4 (bottom panels) has e = 0 and rj = 0.5. The global behavior of S is similar to 
what found for the analytical models of Section 3, with the larger contribution of S near the 
center. In this case we found maximum values of S /~ 0.3. Even larger values were found 
for more battened models: for instance, 0.2 < S/gt^ ^ 0.6 in the central regions (r < Tc) 
of a triaxial model with rj = 0.8, e = 0.6. This is curious because BM95 found similar values 
of S/g^ only in the quite artihcial case of disks with a central hole. This last result suggests 
that when studying MOND equilibrium models for disks or elliptical galaxies it could be 
dangerous to solve equation (7) just relying on the Ansatz that in any case S would be small. 


5. Summary and conclusions 

In this paper we presented a few representative axisymmetric and triaxial density mod¬ 
els, both analytical and numerical, obeying the MOND held equation. The analytical density- 
potential pairs are constructed by means of a simple method based on the deformation of 
the potential of spherically symmetric systems. We show that in this way it is possible to 
build systems with radial prohles and shapes similar to those of real galaxies. Our method, 
although applied in the present paper by using simple deformations is easily generalizable 
to more complicated cases. As a complement to the analytical method, we also presented a 
numerical code (based on spectral methods) to solve the MOND held equation. We tested 
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the code against the new analytical triaxial models and the MOND Kuzmin disk, with ex¬ 
cellent results in terms of both computational time and accuracy of the numerical solution. 
As a hrst simple application of the code, we explored the relevance of the solenoidal field S 
in MOND distributions. Though we conhrm that this held is typically small compared to 
the Newtonian held of the density distribution, we also found that for some systems S is 
certainly not negligible, at least in some regions of space. 

A particularly important application of the numerical solver will be its implementation 
in a particle-mesh N-body code. Such an implementation is promising, because our MOND 
solver is based on an iterative scheme, and at each time in a N-body simulation the potential 
at the previous time is a very good seed for the iterative procedure that should converge 
efficiently. More immediate possible applications of the presented analytical models and 
numerical code are the study of MOND orbits in aspherical density distributions, and the 
vertical motions of stars near the galactic plane. We conclude by pointing out that the 
developed numerical code can also be used to solve the equation for the scalar held in the 
non-relativistic limit of TeVeS, the relativistic MOND theory introduced by Bekenstein (2004; 
equation 53), allowing, for example, to investigate MOND gravitational lensing from non- 
spherical lenses. 

We are grateful to James Binney, Hongsheng Zhao, and the anonymous referee for very 
useful comments. This work was partially supported by a MIUR grant CoFin2004. 


A. Three existence theorems 

We present here three existence theorems that are at the basis of the technique presented 
in Section 2. The hrst result holds in general, while the second and third results assume 
spherical symmetry. 

Theorem 1 Let 0n the Newtonian potential generated by a positive density distribution, 

i.e. 

VVn > 0 (Al) 


over the whole spaee. If 


V-(||V0n||V0n)>O 


(A2) 
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over the whole space, then also 


V- 


/i 


niV0N|h 
\ 0,0 J 


V0N 


> 0 , 


(A3) 


where fi is given in equation (5). 

Proof Equation (A3) can be rewritten as 




Oo 


V0N 


+ S'n 


V 72 , I ■ V0N \ 

+ TT^ j 


(A4) 


from condition (Al) a negative density can be produced only in the region of space where 
Vf/N ■ V0N < 0. If 11“ = 0 the theorem is proved, thus let 11“ ^ 0. In this case, by expansion 
of inequality (A2) we have that on 11“ 


2 , ^ V7 V7A ^ V^-N-V^N 

5'nV 0N > —V5'N-V0N >-- 2 / 2 ’ 

1 + ^N/ao 


(AS) 


and this proves the theorem. 


We now prove a consequence of Theorem 1 holding for spherically symmetric systems, 
which is relevant to the construction of exact aspherical MOND solutions. 

Theorem 2 Let (f){r) the dMOND potential of a spherically symmetric (positive) density 
distribution. Then 


and 


over the whole space. 


/i 


^ (f > 0, 

^^•#>11 


V 


V0 


> 0 , 


(A6) 

(A7) 


Proof Direct calculation of in spherical symmetry, as done in equation (9), proves 
equation (A6). Conditions (Al) and (A2) are then verified and this proves equation (A7). 

Note that the converse of Theorem 2 is not true, i.e., > 0 is not a sufficient condition 

for V ■ (IIV0IIV0) > 0. In fact, the following result holds 


Theorem 3 Let 


Then 


= 47rGp(r) > 0. 
V-(||V0||V0)>O, 


(AS) 

(A9) 
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if and only if 


p{r) > 


M{r) 

47rr3 


Vr. 


(AlO) 


Proof The proof is obtained by expanding inequality (A9) in spherical symmetry as 

df) f d'^cj) 


dr V 


+ AttGp > 0 : 


(All) 


the identity dcf/dr = GM{r)/r‘^ proves the theorem. 
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